Theory of Small Para-Hydrogen Clusters: Magic Numbers and Superfluid Sizes 
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The interplay between magic number stabilities and superfluidity of small para-hydrogen clusters 
with sizes TV = 5 to 40 and temperatures 0.5 K < T < 4.5 K is explored with classical and 
quantum Path Integral Monte Carlo calculations. Clusters with TV < 26 and T < 1.5 K have large 
superfluid fractions even at the stable magic numbers 13, 19, and 23. In larger clusters, superfluidity 
is quenched especially at the magic numbers 23, 26, 29, 32, and 37 while below 1 K, superfluidity 
is recovered for the pairs (27,28), (30,31), and (35,36). For all clusters superfluidity is localized at 
the surface and correlates with long exchange cycles involving loosely bound surface molecules. 



Hydrogen is the simplest and most ubiquitous of all 
molecules in the universe. On earth, it plays an impor- 
tant role in many chemical reactions and is presently be- 
ing developed as an energy transport medium The 
j = rotational state of the para- #2 nuclear spin config- 
uration is, like 4 iJe, a spin-less boson and below about 
6 K has been predicted to be the only naturally occurring 
superfluid besides the helium isotopes 0]. The observa- 
tion of superfluidity in the bulk so far has been thwarted 
by its solidification at 13.96 K. In 1991, Sindzingre et 
al. showed theoretically that small pure clusters with 
13 and 18 molecules were superfluid below about 2 K, 
while larger clusters with 33 atoms had a much smaller 
superfluid fraction Later, Grebenev et al. observed a 
superfluid response in small clusters consisting of 15 — 17 
para Hi molecules surrounding an OCS chromophore all 
within a large helium droplet Q. More recently magic 
cluster sizes of pure pHi clusters with TV — 13, 33 and 55 
were observed with Raman spectroscopy in a cryogenic 
free jet expansion Q . The earlier calculation of superflu- 
idity in pure H2 clusters, as well as the interpretation of 
the OC S{pH2) n experiments have since been confirmed. 
Recently, several groups have reported evidence for magic 
number stabilities at TV — 13@,S[I3| and at 13 and 26 
with a reduced superfluidity [111 ]. 

An intriguing aspect of these studies is the appar- 
ent contradiction between the large superfluid fractions 
and the structured radial distribution functions and the 
magic number stabilities, which indicate a solid-like rigid- 
ity. It is only at high temperatures, when the clusters 
are molten that the radial distributions show the same 
constant interior density and smooth fall-off at the sur- 
face 0, @ found for superfluid helium clusters which are 
known to be liquid [12j. This inconsistency has lead to 
the speculation that pH2 clusters may be considered as 
microscopic supersolids Q. Thus the present study was 
undertaken to clarify how a cluster which appears to be 
solid can also be superfluid. 

To resolve this apparent incompatibility, both classical 
and quantum Path Integral Monte Carlo (PIMC) calcu- 
lations are reported for all sizes between TV = 5 and 40 



for 5 temperatures 0.5, 1.0, 1.5, 3.0 K and the experi- 
mental accessible 4.5 K [l_3j. The calculations indicate 
that clusters with TV = 13, 19, 23, 26, 29, 32, 34 and 37 
have highly symmetric structures and show a propensity 
for stability in agreement with magic number stabilities 
reported earlier for solid ionized rare gas clusters 14 , 15| . 



In the following, these special sizes will be referred to as 
magic clusters. All the H2 clusters with TV < 26 have 
significant superfluid fractions at T < 1.5 K which is 
only slightly suppressed in the magic clusters. At magic 
TV = 26, superfluidity is greatly quenched and for 29, 32, 
34 and 37 reduced even down to T = 0.5 K . Clusters 
with size pairs (24,25), (27,28), (30,31), the singleton 
33, as well as (35, 36) show significantly greater super- 
fluid fractions than their more stable magic neighbours. 
The radial distributions of the superfluid fraction and the 
distribution of permutation cycle lengths reveal that su- 
perfluidity in all the clusters is localized at the surface 
and for the larger superfluid sizes, it correlates with the 
presence of loosely bound surface molecules. These new 
results now clarify the apparent contradiction between 
the structured radial distributions and the large super- 
fluidity found in the previous calculations 0, □J. 

The calculations were carried using the PIMC method 
which is based on the quantum-classical isomorphism 
where each particle is replaced by a polymer made up 
of M "beads" as explained in detail in reference [li} ■ In 
the present quantum calculations a time step of r — 1/80 
was sufficient to obtain converged results within the pair- 
product approximation. Bose statistics is introduced by 
cross linking the polymers to form chains of permut- 
ing cycles (polymers). The classical calculations involve 
no quantum effects such as permutations and are de- 
scribed in Ref. [17|. Two intermolecular potentials were 
used: 1) a Lennard-Jones (LJ) potential with parame- 
ters a — 2. 96 A and e = 34.16 K and 2) the more ac- 
curate Silvera-Goldman (SG) potential. The initial con- 
figurations were chosen either from the Cambridge Clus- 
ter Database [l8| or by carving out a spherical region 
centered around a molecule in an hep H2 crystal. The 
important effects reported here were independent of the 
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potential and the initial configuration. 
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FIG. 1: Energy differences AE as a function of cluster sizes 
calculated classically (a) and with a quantum mechanical 
PIMC program (b) . The sharp minima marked by downward 
arrows are attributed to especially stable magic icosahedral- 
derived structures 114. fl5il . 



Magic classical and quantum clusters were identified by 
examining the energy differences AE = E(N, T) — E(N — 
1,T), where E(N,T) is the total internal energy for a 
cluster of N molecules at temperature T. AE approaches 
the chemical potential at T = 0, hence sizes with AE val- 
ues lower than their neighbours are more localized and 
more tightly bound. Both, the classical and the quantum 
results show about the same enhanced relative stabilities 
at the magic cluster sizes iV = 13, 19, 23, 26, 29, and 32 
indicating that the classical stabilities persist in the quan- 
tum clusters. Although 36 and 38 appear to be "magic" 
in the classical simulation, the expected magic 34 and 37, 
are found only in the quantum calculations (l4, 15 1. All 
our quantum calculated magic numbers agree with well- 
known high symmetry icosahedral-derived structures 14 1 
and have been repeatedly observed in ionized rare gas 
clusters [14j, |l5j. So far, however, only magic N — 13, 33 
and 55, have been observed experimentally Whereas 
N = 13 and 55 correspond to the closing of the first and 
second icosahedral shells, N = 33 is not a magic number 
among the possible modified icosahedral structures but 
has a dodecahedral form. Since the experimental resolu- 
tion H was not sufficient to distinguish from the nearby 
magic 32 and 34, it is possible that the earlier assignment 
may be incorrect. Especially if one considers that dur- 
ing cluster growth icosahedral modified structures will 
be preferred as a result of build-up around the smallest 
magic N = 13 which will be a nucleus for further growth. 

The superfluid fraction p s j Ptotai [IB] is shown in Fig[2] 
as a function of cluster size and temperature. The nearly 
100% fraction in the small clusters decreases sharply 
above 1.5 K |19j |. At each of the magic sizes p s /ptot al- 
most always shows a downward dip reminiscent of those 
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FIG. 2: The total superfluid density fraction is plotted as 
a function of cluster size for five different temperatures. 
The downward arrows indicate the magic clusters with high 
symmetry. The upward pointing arrows indicate islands of 
sizes exhibiting enhanced superfluidity. Superfluidity beyond 
N = 26 above 3 K is below 3%. 



in Fig. 1. These dips are smallest in the small clusters 
but are more pronounced for N > 26. Thus, the classical 
rigidity at the magic sizes suppresses the quantum de- 
localization needed for superfluidity, an effect which be- 
comes stronger with increasing size. Surprisingly, beyond 
magic N = 23, the superfluid fraction at 1 K jumps back 
to about unity for the next larger sizes (24,25). Then 
beyond magic TV = 26 a similar rebound occurs for the 
two lowest temperatures at (27,28). The minimum at 
N = 23 and two maxima at N = 25, 27 were also found 
in recent calculations by Mezzacapo and Bonisegni [ll|. 
In the present calculations similar maxima are found at 
(30, 31), and (35,36). The singleton N = 33 has a very 
weak rebound and seems to follow the behaviour of the 
neighbouring magic sizes as its relatively low binding en- 
ergy in Fig. 1(b) also suggests. Thus in clusters with solid- 
like stable structures superfluidity is suppressed, with the 
largest suppressions found for N > 26, while at sizes not 
corresponding to known magic sizes, it is restored. It is 
interesting to note that at 1.5 K the interactions are so 
strong that superfluidity is almost suppressed for clus- 
ters N > 26, yet quantum derealization in the smaller 
clusters is still sufficient for their superfluidity. 

To understand these unexpected out-of-phase oscillla- 
tions between magic sizes and superfluid sizes the three 
types of radial distributions shown in Fig. 3 were calcu- 
lated for the magic N = 26 cluster and the two adja- 
cent less stable but more superfluid clusters. Fig. 3 (a) 
shows the classical distributions at 1.5 K. As expected 
classicaly(Fig. 3(a)), the molecules in N = 26 follow an or- 
derly partitioning into four localised and distinct groups, 
whereas for N — 25, 27 they are randomly distributed. 
The quantum distributions (Fig. 3b) with only two broad 
peaks show evidence of quantum exchanges between the 
shells. The magic N — 26 cluster has a noticeably larger 
inner peak compared to its neighbors and also shows a 
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FIG. 3: Comparison of radial density distributions for the 
magic cluster N = 26 and the neighbours (25, 27) for three 
characteristic temperatures. The distributions in (a) are clas- 
sical and those in (b) are from quantum calculations. The 
radial distribution of the superfluid fraction is plotted in (c). 
We show one temperature in (a) since the classical distribu- 
tions are less sensitive to temperature. 



slight dimple in the second outermost peak reminiscent 
of the classical distribution. The radial dependence of 
the superfluid fraction p s (r)/p [20j in Fig. 3(c) is com- 
puted by binning the radial location of the beads that 
are involved in permutation cycles [2oT |. These distribu- 
tions exhibit large differences with temperature and with 
sizes as expected from Fig. 2. Surprisingly, however, su- 
perfluidity is not greatest near the center as found for 
4 He clusters 12j, but is localized at the surface beyond 
the outer maximum in the radial density distribution. 
The superfluidity is small in the inner shell specially for 
magic N = 26. We note that the apparent randomness in 
the classical radial curves (Fig. 3(a)), which is suggestive 
of less rigid and symmetric structures for the non-magic 
N = 25 and 27, correlates with the increased superfluid 
fractions (Fig. 3(c) and Fig. 2). 

Additional insight comes from Fig. 4 which provides a 
cross section view of N = 13, 26 and 27 and the probabil- 
ity distributions of the permutation cycles as a function 
of the permutation lengths. The contours for N = 13 
show a considerable delocalization of the outer layer. The 
A^ = 13 permutation probabilities are dominated by cy- 
cles with lengths of 5 and 6 molecules, corresponding to 
rings around the center, with a relatively high probability 
of 10% in accordance with the large superfluid fraction 
(Fig. 2). Even though the central molecule appears to be 
localized in the contour plot at 0.5 K, it also participates 
in ring exchanges as indicated by the small « 0.1% prob- 
ability for cycle exchange lengths of 13. The contour 
plot for magic N = 26 indicates that all its molecules 
are highly localized and that it is solid-like justifying the 
greatly suppressed superfluidity (Fig. 2). Contours at 
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FIG. 4: Top panels: Density contour plots at z w plane for 
magic N = 13, 26 and non-magic N = 27. Bottom panels:The 
probability distributions for permutations of a given length. 
The contours for magic N — 13, show delocalized molecules 
at the outer layer that take part in cycles with lengths 5 and 
6 which correspond to rings around the central molecule. In 
magic N — 26 at 1.0 K, all the molecules are highly localized. 
The surface molecules in N — 27 cluster at 1 K are highly 
delocalized, whereas the core molecules are nearly rigid. This 
correlates with the long exchange paths around the localized 
core shown below the contour plot. 



different values of z for larger clusters at T > 1.5 K (not 
shown) show similar localization. The corresponding per- 
mutation probabilities are peaked at very low cycle num- 
bers and are mostly less than 0.05%, which explains their 
small residual superfluidity in Fig. 2. But as shown by 
the permutation probabilities even the magic A^ = 26 
melts and become superfluid at 0.5 K with a predomi- 
nance of cycles from 7 to 15. The top panel in Fig. 4 
for non-magic N = 27 has a liquid outer layer similar to 
magic 13 in agreement with Fig. 3(c) but with a smaller 
peak permutation probability of about 3.5% as expected 
from its small overall superfluid fraction (Fig. 2). At 
0.5 K the permutation probability curve smoothens and 
extends out to include cycles equal to 27 suggesting that 
also the core molecules are participating in the permu- 
tations. Thus both the contour plots and probability 
distributions of the N — 21 cluster are consistent with 
the onset of superfluidity in the surface region. 

Recently, Mezzacapo and Boninsegni also observed an 
enhanced superfluidity for N — 25, but their conclusion 
that "the addition of a molecule to the A^ + 1 = 26 has 
the effect of frustrating the solid order of the inner shell, 
increasing molecule delocalization and leading to quan- 
tum exchanges" ll| is at variance with these new results, 
which clearly show that superfluidity is at the surface. 

The nature of the disorder favouring a large superfluid 
fraction at the surface emerges from the "inherent struc- 
ture" (IS) analysis of Stillinger and Weber [21]] : "inherent 
structures which underlie the liquid state are those sta- 
ble particle packings (potential minima) which can be 
reached by a steepest descent quench on the potential 
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energy hypersurface" . This quenching procedure elimi- 
nates all kinetic effects due to thermal excitations or zero- 
point motion. To generalize this concept to quantum me- 
chanical systems, in the steepest descent minimization, 
the gradient was calculated from the path integral action 
rather than from the potential. The permutations were 
turned off and the IS analysis was applied to "Boltzman- 
nons" . The quenched structures are, in general, indepen- 
dent of the temperature 2l| and are identical to Wales 
classical minimum potential energy configurations [lil ]. 
This geometric spatial correspondence explains the ori- 
gin of the persistence of the classical structures in the 
energy and superfluid densities discussed above. 

Among the small cluster sizes, a quenched configura- 
tion similar to the Wales body centered icosahedron was 
quite often found for N — 13 indicating that it is par- 
ticularly stable. Occasionally, quenching would generate 
variant structures close to the classical ones with some 
delocalized molecules on the surface suggestive of melt- 
ing. Clusters differing by one or two molecules from the 
magic sizes are more often seen to have defective sur- 
faces. Since these clusters differ essentially only in the 
bonding of the outer molecules, their smaller binding en- 
ergies (Fig. 2) indicate that these outer molecules are less 
tightly bound than in the case of the magic clusters. For 
example, N — 18 would statistically appear more often 
with structures that deviate slightly from the classical 
clusters, while for magic N = 19, almost every IS cluster 
is the same as the classical cluster. The IS analysis for the 
other larger superfluid sizes indicates that their surface 
molecules are also less tightly bound and less ordered. 

In summary, our analysis reveals that pure pH2 clus- 
ters with N < 26 at temperatures T < 1.5 K are liquid- 
like and have a large superfluid response. It is only some- 
what reduced in magic clusters N = 13, 19 and 23, which 
are classical magic sizes with highly symmetric icosahe- 
dral structures. According to Fig 1 (b), the difference in 
internal energies, which is the energy needed to add one 
molecule, is less than 38 K for these highly superfluid 
clusters. The larger magic clusters N = 26,29,32,34, 
and 37, in which superfluidity is strongly quenched at 
temperatures T > 0.5 K, all have considerably larger 
internal energy differences of more than about 40 K . 
Superfluidity is restored in the cluster size pairs N — 
(24,25), (27,28), (30,31), the singleton 33, and (35,36), 
with smaller internal energy differences compared to the 
magic clusters which lie in between. In these superfluid 
sizes, quantum derealization of the loosely bound ad- 
molecules enables them to explore many different surface 
structures, thereby favouring large permutation cycles 
and an increased superfluidity. Our calculations reveal 
that with increasing cluster size the strong many-body 
intermolecular interactions lead to a rigid solid-like inner 
core thereby pushing the derealization induced super- 
fluidity towards the surface, where it is favored by the 



reduced coordination and weak inward interaction with 
the small central core. The overall decay of superfluid- 
ity with cluster size and its increased localization on the 
surface agrees with the macroscopic limit of 2D surface 
superfluidity and zero response in the bulk [22l |. 
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